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Abstract. Static and dynamic properties of two-dimensional bidisperse dissipative particles are numerically studied near the 
jamming transition. We investigate the dependency of the critical scaling on the ratio of the different diameters and find a 
new scaling of the maximum overlap (not consistent with the scaling of the mean overlap). The ratio of kinetic and potential 
(J_) ' energies drastically slows down near the jamming transition, i.e. the relaxation time diverges at the jamming point. 
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INTRODUCTION 

Jamming is an universal feature of both thermal and athermal systems, for instance, glasses, granular particles, 
emulsions, colloidal suspensions, foams, etc. where constituents are arrested in disordered states so that the material 
obtains rigidity. The jamming transition is governed by temperature, density and external loads, and a lot of systems 
can be mapped onto a unified phase diagram [1,2, 3]. 

Jamming of athermal systems, i.e. granular particles [4, 5, 6, 7], emulsions [8, 9] and foams [10, 1 1], occurs at zero 
temperature at the critical density (area fraction) <j> c [1]. At this point, each particle begins to touch its neighbors and the 
mean coordination number z, defined as the average number of contacts per particle, jumps from zero to the isostatic 
value z c = 2d in if -dimensions [12]. Some macroscopic variables also indicate the acquisition of rigidity, for example, 
the pressure p and the shear modulus G start to increase from zero, and the bulk modulus K discontinuously jumps to a 
non-zero value [4, 5, 6, 7]. Above this threshold, the excess coordination number Az = z — Z c , p, G and K respectively 
scale as Az ~ A<j) '/ 2 , p ~ A<j) V, G ~ A0 y and K ~ A0^ with the distance from the jamming point A<j) = <j> — (j) c . In the 
case of monodisperse particles, the first peak of the radial distribution function gi scales as g\ ~ A0 [13, 14, 15, 16]. 
In the case of bidisperse particles, one can also see a similar divergence of the radial distribution function with the 
scaled distance [7]. Moreover, the mean overlap between particles (8) linearly scales as (8) ~ A0. 

Those power law scalings are analogous to those found in critical phenomena. However, some variables show 
discontinuous changes at the critical point as in first-order phase transitions. Furthermore, the critical exponents y/, y 
and X depend on the interparticle forces, which suggests that the jamming transition is entirely different from usual 
critical phenomena [4, 5, 6, 7]. 

Even though the critical scalings above the jamming transition are well established, not much attention has been 
paid to the critical amplitudes. Many previous works on bidisperse systems only focused on the particular case that the 
ratio of the different diameters equals p = 1 .4. Furthermore, the dynamic properties are not understood yet. 

In this paper, we study the static and dynamic properties of two-dimensional bidisperse particle systems by 
numerical simulations and investigate the dependency of the critical amplitudes on p. We also study the dynamic 
properties near the jamming transition, where we adopt the ratio of kinetic to potential energies to quantify relaxation 
and show its drastic slowing down near the jamming point. 



MOLECULAR DYNAMICS SIMULATION 

We investigate two-dimensional packings of bidisperse dissipative particles using molecular dynamics (MD) simula- 
tions. Our strategy is to increase the diameter of each particle until a desired area fraction <p is obtained. After the 
desired area fraction is obtained, the kinetic energy decreases to zero due to the dissipative forces between particles 
and the system relaxes to its static state. We study dynamic properties during the relaxation and static properties after 
the relaxation. 



At first, we prepare a binary mixture of large and small particles with initial diameters <7l(0) and Os(0), respectively, 
where the total number of particles is N = 32768. We randomly distribute them in a L x L periodic box one-by-one, 
avoiding overlap between particles. Then, we slowly increase the diameter of each particle [17, 18, 19, 20] as: 

a a (t)=g r a a (0) (a = LorS) (1) 

with a constant growth rate g r . Because we fix the mass density of particles, the mass also increases while we 
increase the diameter O a (t) = G a (0)(g r t + 1). Therefore, the size ratio p = <7l(0/ c7 s(0 = Ol(0)/os(0) does not 
change throughout our simulations. 

Each particle ; can be in contact with several other particles j and moves according to the equation of motion 

mi {t)ti = Y J { k " S 'J _ 77,1 ( V <J ' n 'i ) } n 'J ' ^ 

j 

where m,(f), k n , r\ n , mj = (x; — Xy)/|Xj — X;|, and V(j = x,- — Xj are the mass of particle i, the spring constant, the 
viscosity coefficient, the overlap, the normal unit vector, and the relative velocity, respectively. The vector X; is the 
position of particle i, and each dot above x represents the time derivative. 

When the area fraction reaches the desired value (j) at t = to, we stop to increase the diameter of each particle. Then, 
the kinetic energy of the system decreases due to the dissipative forces and we assume the system is relaxed to its static 
state if the ratio of kinetic to potential energies becomes lower than 10~ 6 . In the following, we scale mass, length and 
time by m u = m^ito), l„ = C7l(?o) and f„ = m^to) / rj n , respectively, and use g r = 10~ 4 /f„, k n — 4.0 x 10 5 m H /f 2 and 
r\n = m u /t u . 

From Eq. (2), the restitution coefficient in the relaxation stage (t > to) defined as the ratio of speeds after and before 
a collinear collision is calculated as [21] 

e„ = e-™° , (3) 

where t c = %j J (k n /mij) — ?7q and Tjo = T\ n / (2m/y) are a typical response time and the rescaled viscosity coefficient 

with the reduced mass mij = m,'(fo)my(fo)/(m,(fo) +mj(to)). Although the restitution coefficients between two small 
particles, two large particles, and small and large particles are slightly different, we find e n = 0.99 ±0.003 in the three 
cases since dissipation is rather weak. 



CRITICAL SCALING 

Changing the values of <p and p between 0.8 < (f) < 0.9 and 1 .2 < p < 2.4, respectively, we repeat our simulations and 
measure the mean overlap (8), the mean coordination number z, the pressure p, the first peak of the radial distribution 
function with the scaled distance g + and the maximum overlap 8 m , after each system is relaxed to its static state. We 
define the critical area fraction <f> c as the point at which (8) starts to increase from zero and z jumps from zero to the 
isostatic value z c — 4. It should be noted that we remove rattlers of which the number of contacts are less than 3 from 
the statistics, because they do not contribute to the contact- and force-chain networks [4, 5, 6, 7]. 

The critical area fraction depends on both p and g r . Figure 1 shows (j) c as a function of the relative standard deviation 

R = (a) ' (4) 

where (a) = (1 +p)(7s(fo)/2 and (a 2 ) = (1 + p 2 )as(fo) 2 /2 are the mean diameter and the mean square of diameter, 
respectively. In this figure, the closed circles are our results and the open circles are the results of another simulation 
of bidisperse hard spheres in two dimensions [22], where a large discrepancy between these results can be seen below 
R = 0.2, mainly caused by the dependency of <j) c on g r . In this paper, we do not discuss the rate dependency of <j) c , 
however, it is known that <j) c strongly depends on g r , if R is small (see Ref. [23] for a more detailed discussion). We 
also list our results of (j) c for different p and R in Table 1 . 
Above (j> c , (8) scales as [4, 5, 6, 7] 

<S)=B 5 (p)Af, (5) 

where pi and B§(p) are the critical exponent and the critical amplitude, respectively (see Table 1). From our simula- 
tions, we find jj. ~ 1.0, i.e. (8) depends linearly on A0 (see Table 2). 



o 

0.88 - 

(f>c 

o 

0.86 - 0« 

o # o • 

°%oq °* 

0.84 I 1 1 

0.2 0.4 0.6 

R 

FIGURE 1. Critical area fraction <j) c as a function of R, where the closed circles are the results of our simulations and the open 
circles are the results of another simulation of bidisperse hard spheres in two dimensions [22]. 



TABLE 1. The critical area fraction and the critical amplitudes. 



p 


R 




B 5 






A+ 


Am 


A m /B s 


1.2 


0.127 


0.846054 


0.262691 


0.062944 


9.623507 


0.108473 


3.271311 


12.453076 


1.4 


0.230 


0.847454 


0.247794 


0.062175 


9.789582 


0.112457 


3.839339 


15.494075 


1.6 


0.313 


0.851053 


0.234816 


0.058778 


9.696656 


0.107054 


4.508161 


19.198696 


1.8 


0.381 


0.853553 


0.240875 


0.061023 


10.680258 


0.102261 


4.554234 


18.907043 


2.0 


0.435 


0.855449 


0.246302 


0.062333 


9.661944 


0.109497 


4.378124 


17.775430 


2.2 


0.481 


0.858054 


0.278428 


0.067636 


9.986629 


0.122986 


5.144424 


18.476676 


2.4 


0.518 


0.859950 


0.268369 


0.073998 


10.649745 


0.207261 


5.304872 


19.767081 



We define the pressure as the virial pressure p = Y,i<j r ijfij/L 2 with the interparticle distance r,,- and force fy [24] 
and introduce the excess coordination number as Az = z — z c - Because we use bidisperse particles, the usual radial 
distribution function has three peaks around as, 0L and (cl + Os)/2, respectively. However, if we introduce the scaled 
distance between particles i and j as £, = r, 7 / (r, + rj), where r, and rj are the radii of particles i and j, respectively, 
the radial distribution function of | has the first peak g + around % = \. Then, we find p, Az and g + scale as 

P = A p (p)(5)V, (6) 

Az = A z (p)(8)^ (7) 
g + = A + (p)(8)-^, (8) 

respectively, where A p (p),A z (p) and A + (p) are the critical amplitudes (see Table 1), and i/^~l,£~l/2 and 77+ ~ 1 
are the critical exponents (see Table 2). Figures 2 (a)-(c) display p* = p/A p (p), Az* = Az/A z (p) and g* + = g+/A+(p), 
respectively. Since (8) ~ A0, these results confirm p ~ A0, Az ~ A0 1 / 2 , g+ ~ A<j)~ l and g+(8) w const [4, 5, 6, 7]. 
From our simulations, we also find 5m scales as 

5n=A m (p)(5) A (9) 

with the critical amplitude A m (p) (see Table 1). Figure 2 (d) shows 5jjj = <5 m /A m (p), where we find X ~ 1 similar to 
(see Table 2). The ratios A m (p)/Bg(p) are not constant (see Table 1), which means that the probability distribution 
functions of the overlap are not self-similar, but change shape with increasing p [25]. 
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FIGURE 2. (Color online) Double logarithmic plots of (a) p*, (b) z* , (c) g\ and (d) 5,* as functions of (8). The solid circles, 
open circles, solid squares, open squares, crosses, open triangles and solid triangles are the data for p = 1.2, 1.4, 1.6, 1.8, 2.0, 2.2 
and 2.4, respectively. 



TABLE 2. Estimated exponents and errors. 



function 


base 


exponent 


value 


deviation 


error [%] 


(*) 


A<j> 




1.006 


±0.003 


0.349 


P 


(8) 


¥ 


1.039 


±0.001 


0.065 


Az 


(5) 


c 


0.551 


±0.001 


0.124 


g+ 


(S) 


'?+ 


0.960 


±0.001 


0.148 


Sm 


(8) 




0.976 


±0.001 


0.164 



CRITICAL SLOWING DOWN NEAR JAMMING 



After we stop to increase the diameter of each particle at t = to, the system relaxes to its static state. To quantify the 
relaxation dynamics, we introduce a dimensionless energy 



X(t) 



K(t Q )/U(t ) ' 



(10) 



where K(t) and U(t) are the kinetic and potential energies, respectively. The ratio K(t)/U(t) can be used to estimate 
<p c , because U(t) drops to zero and K(t)/U(t) diverges if < (j> c [26]. The function %(t) is the ratio K(t)/U (t) scaled 
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FIGURE 3. (Color online) (a) %(t) against logarithmic time above <j) c , where the circles, triangles, squares and diamonds are the 
results of simulations with A<j> = 1.2 x 10~ 2 , 2.5 x 10~ 3 , 4.4 x 10~ 4 and 1.4 x 1CT 4 , respectively. The dotted lines represent Eq. 
(11). (b) Double logarithmic plot of %(t)/Ci, where the circles, triangles, squares and diamonds are the results of simulations with 
A(j) = 5.5 x 10~ 4 , 4.5 x 10~ 4 , 2.5 x 10~ 4 and 1.5 x 10~ 4 , respectively. The dotted line represents the power law decay with the 
exponent 2.541. Here, we used p = 1.4. 
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FIGURE 4. (Color online) (a) T and (b) |3 against logarithmic (8), where the insets show the double logarithmic plots of them. 
The solid circles are our simulation results and the solid lines represent T °c {8}~ K and /3 oc (5}~ v with K = 0.274 ±0.006 and 
v = 0.1 12 ±0.008, respectively. Here, we used p = 1.4. 



by the value at ?o. In the following, we redefine to — for simplicity and study the case of p = 1 .4. 

Figure 3(a) shows our results of x(t) against logarithmic time, where the decrease of %(t ) drastically slows down 
near the jamming point. To quantify the slowing down, we fit %(t) by a stretched exponential function 

X(rH^ ( ' /T)/ \ (ID 

where T and j3 are the relaxation time and the stretching exponent, respectively. Figure 4 displays both quantities as 
functions of (5). Slightly above e (0 < A0 < 1CT 3 ), both T and j3 scale as T ^ (8)~ K and j3 °< (<5}~ v with K ~ 0.27 
and v ~ 0.1 1, respectively. However, given the rather narrow range of T and /3, power laws with such small exponents 
are not of significance and require a wider range before conclusions can be drawn. It should be noted that j3 > 1 near 
the jamming point and %{t ) decays faster than exponential with time. 

Figure 3(b) is a double logarithmic plot of the long term asymptotic behavior of % (t ) near the jamming point, where 
we find a power law decay 

x(t) = c L r a . (12) 



The exponent a ~ 2.5 is not significantly changed by changing (j) and faster than the decay rate of the kinetic energy 
in the homogeneous cooling state of granular gases, where the exponent is given by 2 [27]. 

SUMMARY 

In summary, we numerically investigated the static and dynamic properties of two-dimensional bidisperse particles 
near the jamming transition and systematically studied the dependency of the critical exponents and amplitudes 
on the size ratio p. For different size ratios, we found different scaling prefactors of the average and maximum 
overlaps confirming their different behaviors and thus indicating different shapes of the overlap probability distribution 
functions [25]. 

Concerning the dynamics, we used the energy ratio %(t) to quantify the relaxation of kinetic energy and report that 
it resembles a stretched exponential and drastically slows down near the jamming point. The asymptotic behavior of 
X{t ) resembles a power law decay with exponent 2.5, which is faster than the decay rate of the kinetic energy in the 
homogeneous cooling state of granular gases. 

The cases of monodisperse and polydisperse particles are left to future study, as is the case for three dimensions. 
Although the number of particles in our simulations may be large enough to study the critical area fraction and the 
critical scaling, we have to investigate the influence of the system size, the growth rate and the restitution coefficient 
elsewhere. 
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